function [Q,theta2,params] = induceStatVAR_threshold_macro(paramsValues,input)

% We turn paramsValuesinto a struct 
for i=1:size(input.selectParams,1)
    name = input.selectParams(i,1);
    paramsInput.(name{1}) = paramsValues(i,1);
end

% Accounting for calibrated coefficients, we unfold paramsInput by to get params
params = paramsInput;
% The setdiff implies that we do not overwrite elements in params based on calibrateparams
if ~isempty(input.calibrateParams)
    namesCalibrate = setdiff(fieldnames(input.calibrateParams),input.selectParams);
    for i=1:size(namesCalibrate,1)
        name = namesCalibrate(i);
        params.(name{1}) = input.calibrateParams.(name{1});
    end
end

% We unfold theta2
nx                   = size(input.varXData,1);
[~,~,hx_1,hx_2,~,gama0,gamaz,gamax,sigmaz] = unfoldTheta2_threshold_macro(input.theta2,nx);
% Scaling hx_1 and hx_2
%hx_1                = params.delta*hx_1;
%hx_2                = params.delta*hx_2;
hx_1                = params.delta1*hx_1;
hx_2                = params.delta2*hx_2;
if input.hxRegimeOn == 0
    hx_2            = hx_1;
end

% Re-estimating h0_1, h0_2 and sigma
T        = size(input.econAct,1);
regime_1 = input.econAct >= input.thresholdLevel; %a boom
regime_2 = input.econAct < input.thresholdLevel;  %a recression
if input.h0RegimeOn == 1
    h0_1 = (eye(nx)-hx_1)*sum(input.factors(1:nx,:).*repmat(regime_1',nx,1),2)/sum(regime_1);
    h0_2 = (eye(nx)-hx_2)*sum(input.factors(1:nx,:).*repmat(regime_2',nx,1),2)/sum(regime_2);
else
    h0_1 = mean(input.factors(1:nx,:),2)-hx_1*sum(input.factors(1:nx,:).*repmat(regime_1',nx,1),2)/T ...
                                        -hx_2*sum(input.factors(1:nx,:).*repmat(regime_2',nx,1),2)/T; 
    h0_2 = h0_1;
end
Wxhat               = input.factors(1:nx,2:T) - (repmat(h0_1,1,T-1) + hx_1*input.factors(1:nx,1:T-1)).*repmat(regime_1(1:T-1)',nx,1) ...
                                              - (repmat(h0_2,1,T-1) + hx_2*input.factors(1:nx,1:T-1)).*repmat(regime_2(1:T-1)',nx,1);
VarWxHat             = zeros(nx,nx);
for t=1:T-1
    omegaHat = input.VarU(1:nx,1:nx,t+1)...
        + hx_1*input.VarU(1:nx,1:nx,t)*regime_1(t,1)*hx_1' ...
        + hx_2*input.VarU(1:nx,1:nx,t)*regime_2(t,1)*hx_2' ...
        - input.Cov10U(1:nx,1:nx,t+1)*regime_1(t,1)*hx_1' - hx_1*input.Cov01U(1:nx,1:nx,t+1)*regime_1(t,1) ...
        - input.Cov10U(1:nx,1:nx,t+1)*regime_2(t,1)*hx_2' - hx_2*input.Cov01U(1:nx,1:nx,t+1)*regime_2(t,1);    
    VarWxHat = VarWxHat + 1/(T-1-2*(nx+1))*Wxhat(:,t)*Wxhat(:,t)' - 1/(T-1)*omegaHat;
end
[sigma,checkWx] = chol(VarWxHat,'lower');
if checkWx ~= 0
    %sigma = sigmaOld;
    VarWxHat             = zeros(nx,nx);
    for t=1:T-1
        VarWxHat = VarWxHat + 1/(T-1-2*(nx+1))*Wxhat(:,t)*Wxhat(:,t)';
    end
    sigma = chol(VarWxHat,'lower');
end
                                       
eigValues_1         = eig(hx_1);
modulus_1           = sqrt(real(eigValues_1).^2 + imag(eigValues_1).^2);
eigValues_2         = eig(hx_2);
modulus_2           = sqrt(real(eigValues_2).^2 + imag(eigValues_2).^2);
if any(modulus_1>=1) || any(modulus_2>=1)
    Q = 1D35;
else
    % We have a stationary model
    if input.hxRegimeOn == 0 && input.h0RegimeOn == 0
        % Unconditional variance is available in closed-form
        h0_1  = (eye(nx)-hx_1)*mean(input.factors(1:nx,:),2);
        T      = size(input.factors,2);
        residualsTmp = zeros(nx,T);
        for t=2:T
            residualsTmp(:,t) =  input.factors(1:nx,t) - (h0_1 + hx_1*input.factors(1:nx,t-1));
        end
        varU     = sum(input.VarU(1:nx,1:nx,1:T-1),3);
        cov10U   = sum(input.Cov10U(1:nx,1:nx,2:T),3);
        cov01U   = sum(input.Cov01U(1:nx,1:nx,2:T),3);
        omegaTmp = 1/(T-1-nx-1)*(residualsTmp*residualsTmp') + 1/(T-1)*(-varU - hx_1*varU*hx_1' + cov10U*hx_1' + hx_1*cov01U);
        [sigmaTmp,check] = chol(omegaTmp,'lower');
        deScaling     = 1;
        while check ~= 0
            deScaling = deScaling*0.9;
            omegaTmp = 1/(T-1-nx-1)*(residualsTmp*residualsTmp') + 1/(T-1)*deScaling*(-varU - hx_1*varU*hx_1' + cov10U*hx_1' + hx_1*cov01U);
            [sigmaTmp,check] = chol(omegaTmp,'lower');
        end
        varX     = reshape((eye(nx^2)-kron3(hx_1,hx_1))\reshape(sigmaTmp*sigmaTmp',nx^2,1),nx,nx);
    else
        % The unconditional variance is computed by simulation
        numSim         = input.numSim;
        thresholdLevel = input.thresholdLevel;
        rng(1,'twister')
        shocks         = randn(nx,numSim);
        shocksz        = randn(1,numSim);
        xSim           = zeros(nx,numSim);
        zSim           = zeros(1,numSim);
        for s=2:numSim
            if zSim(1,s-1) >= thresholdLevel
                xSim(:,s) = h0_1 + hx_1*xSim(:,s-1) + sigma*shocks(:,s);
            elseif zSim(1,s-1) < thresholdLevel
                xSim(:,s) = h0_2 + hx_2*xSim(:,s-1) + sigma*shocks(:,s);
            end
            zSim(1,s) = gama0 + gamaz*zSim(1,s-1) + gamax'*xSim(:,s-1)+ sigmaz*shocksz(1,s);
        end
        varX = zeros(nx,nx);
        for i=1:nx
            varX(i,i) = var(xSim(i,input.burnin+1:numSim),1,2);
        end
    end
    % Matching all the variances
    Q = 0;
    for i=1:nx
        Q = Q + ((varX(i,i)-input.varXData(i,i))/input.varXData(i,i))^2;
    end
    Q = sqrt(Q/nx)*100;
end

% We save the potentially scaled elements in theta2
theta2 = setTheta2_threshold_macro(h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz);

end

